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O ABSTRACT 

■ We suggest an explanation for the twin kilohertz quasi-periodic oscillations 

(kHz QPOs) in low-mass X-ray binaries (LMXBs) based on magnetohydrody- 
namics (MHD) oscillation modes in neutron star magnetospheres. Including 
the effect of the neutron star spin, we derive several MHD wave modes by 



> 

■ solving the dispersion equations, and propose that the coupling of the two 



resonant MHD modes may lead to the twin kHz QPOs. This model naturally 
relates the upper, lower kHz QPO frequencies with the spin frequencies of the 



neutron stars, and can well account for the measured data of six LMXBs. 

O 
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1 INTRODUCTION 

The fastest variability components in X-ray binaries, the kilohertz quasi-periodic oscillations 
(kHz QPOs), have been detected in about thirty neutron star low-mass X-ray binaries (NS 
LMXBs) since first discovered with Rossi X-ray Timing Explorer in 1996 (see van der Klis 
2006 for a review). Twin kHz QPOs appeared simultaneously in about twenty NS LMXBs. 
These QPOs are thought to reflect the motion of matter at the inner edge of an accretion 
disc around the NS. Early observations showed that the frequency separation (Au) of the 
twin kHz QPOs is close to the NS spin frequency (u s ) (e.g. Strohmayer et al. 1996; Ford et al. 
1997), suggesting a beat-frequency explanation (Miller, Lamb & Psaltis 1998). However, the 
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more detailed measurements revealed that Av is generally inconsistent with a constant value 
of v s but varying with the upper {y-i) or lower iv\) frequency of the twin kHz QPOs (van 
der Klis et al. 1997; Mendez et al. 1998, 1999). Stella & Vietri (1999) propose a relativistic 
precession model which predicts a changing peak separation of the twin kHz QPOs. The 
upper and the lower kHz QPO frequencies are identified with the Keplerian frequency of the 
rotational plasma flow at the inner edge of the disc and the periastron precession frequency 
of the orbit, respectively. In this model, the NS spin frequency plays no role in setting up 
the frequencies of the kHz QPOs, but a massive (~ 2M ) NS is usually required to match 
the observations. Abramowicz et al. (2001, 2003) suggest that the twin kHz QPOs can be 
explained by a nonlinear resonance in the epicyclic motion in the accretion disc, leading to 
the 3 : 2 ratio of the upper and lower frequencies, although whether there is an intrinsically 
preferred ratio between the v 2 and v x is controversial (Belloni et al. 2007). 

The subtle feature of twin kHz QPOs is that Av is variable but seems to be around 
either v s or v s /2 (Smith, Morgan & Bradt 1997; Wijnands & van der Klis 1997; Markwardt, 
Strohmayer & Swank 1999; Chakrabarty et al 2003; Wijnands et al. 2003; Barret, Boutelier & 
Miller 2008; see, however, Yin et al. 2007; Mendez & Belloni 2007). This has motivated some 
scenarios of kHz QPOs taking into account the effect of the NS spin. Osherovich & Titarchuk 
(1999) suggest that kHz QPOs could be modeled as oscillations of the blobs thrown into 
magnetosphere from the inner edge of the accretion disc. The lower kHz QPO frequency v\ 
is identified as the Keplerian frequency in the disc and the upper frequency v 2 the hybrid 
frequency about the Keplerian frequency and the NS spin frequency. Lee, Abramowicz & 
Kluhiak (2004) show that the kHz QPOs may be attributed to forcing of epicyclic motions 
in the accretion disc by the NS, which induces resonance at selected frequencies when the 
frequency separation Av is equal to v s or u s /2. Li & Zhang (2005, see also Zhang 2004) 
present an alternative interpretation for the origin of the twin kHz QPOs by considering the 
interaction between the NS magnetic field and the surrounding accretion disc, which gives 
rise to MHD loop oscillations at the inner edge of the disc with the frequencies depending 
on the spin frequencies. 

In this paper, we propose a model to explain kHz QPOs in NS LMXBs based on the 
interaction of accreting plasma with the NS magnetosphere. This model is partially based 
on the interpretation of Zhang (2004) and Rezania & Samson (2005). In the latter work it 
was argued that distortion of the NS magnetosphere by the infalling plasma of the Keplerian 
accretion flow can excite resonant shear Alfven waves in a region of enhanced density gra- 
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clients, where accretion material flows along the magnetic field lines in the magnetosphere. 
The major difference between our model and Rezania & Samson's (2005) is that we have 
included the effect of the gravity and spin of the NS on the field line resonances. 

This paper is organized as follows. In section 2 we introduce the basic physical model and 
derive the relation between u s , V\ and z^. In section 3 we compare the theoretical relations 
with the observational data of six sources (4U 0614+09, 4U 1608-52, 4U 1636-53, 4U 
1728—34, 4U 1915—05, and XTE 1807—294). In section 4 we summarized our results and 
discuss their possible implications. 



2 THE MODEL 

In LMXBs the plasma from the donor star accretes onto the NS via an accretion disc. Mate- 
rial in the disc firstly rotates in a Keplerian motion, then corotates with the magnetosphere 
after it is trapped completely by the NS magnetic field at the magnetospheric radius, and 
finally flows along the field lines to the polar cap. Some resonant modes may be excited by 
the perturbations at the magnetospheric radius when the plasma begins to corotate with 
the magnetosphere (Osherovich & Titarchuk 1999; Lee et al. 2004; Zhang 2004; Rezania & 
Samson 2005). We consider the QPOs as a modulation effect of the MHD waves which are 
produced at the magnetospheric radius, and the coupling of the two resonant MHD modes 
may lead to the twin kHz QPOs in the power spectrum. 

We consider the MHD equations in a frame of reference corotating with the NS (shown 
in Fig. 1), written as follows (Landau & Lifshitz 1976) 

p — = -VP + J x B + 2pv x ft + pVt x (r x Q) - p^r, (1) 

(J. L 

— = V x (v x B) = (B ■ V)v - (v ■ V)B - (V • v)B, (2) 

f t +V-(pv) = 0, (3) 

Pp~~' = const, (4) 

where v is the plasma velocity, J electric current, B magnetic field, r the displacement from 
the center of the NS to the plasma, p plasma density, P barometric pressure, 7 adiabatic 
index, G gravitational constant, M and ft the mass and the angular velocity of the NS^, 

1 Actually ii is the angular velocity of the NS magnetosphere, which could be slightly deviate from that of the NS. 
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respectively. The third, fourth and fifth terms on the rhs of Eq. (1) represent the Coriolis 
force, the centrifugal force and the gravity, respectively. 

Observationally the accretion rates in LMXBs change on a timescale of ~ 10 3 — 10 4 s. 
This is much more than the relaxation (or dynamic) timescale of the plasma at the inner 
disc radius (~ 1CT 2 — lCT 3 s). So we can approximate the plasma to be in an equilibrium 
state, which is subject to small perturbations, as discussed in Benz (2002). By use of the 
current expression, 

J = x B, (5) 

where \i is is vacuum magnetic conductivity, Eqs. (l)-(3) can be transformed to be 

Po^ + («o-V)u = -VP + ±(VxB o )xB o + 2p v o xn+p n 2 r o -p (n-r o )n-p ^fr o ,(6) 
dB 



dt 
dpo 



[B ■ V)v - (v ■ V)B - (V • Vo)B , (7) 

m + V • (pv ) = 0, (8) 

where the subscript denotes variables in the equilibrium state. The initial relative velocity 
(vo) is equal to zero in the corotating reference system and can be expressed as O x r in 
the inertial reference system. 

Now we consider the MHD equations for the plasma subject to small perturbations, 

p f + (v ■ V)v = -VP + i(V x B) x B + 2p G v x ft + p n 2 r - p (n ■ r)Q - p ^r, (9) 
f)B 

—— = (B ■ V)v — (v ■ V)B — (V • v)B. (10) 

f f +V- (pv) = 0, (11) 

Pp-^ = P oPo \ (12) 

where v = v + v s = v s , B = B + B s , r = r + r 8 , p = p + p s , P = P + P s with the 
subscript s denoting the perturbed quantities (v s <g; \ fl x r \, B s <^ B , r s <^ r , p s <^ p , 
P s P ) and with the superscript " the variables after the disturbance. Combining Eqs. (6)- 
(12) we get the equations about the perturbed quantities in the first order approximation, 

Po^f = -^VPs + l[(Vy<B )xB s + (VxB s )xB ]+2p v s xn + p n 2 r s 
-p (n.r s )n-p ^-r s , 

= (B ■ VH - (V • v s )B , (14) 



and 
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^± + p V-v s = 0. (15) 

Differentiating Eq. (13) and substituting V(B ■ B s ) = (B ■ V)B S + (B s ■ V)B + B s x 
(V x B ) + B x (V x B s ) into it give 

= -^|(vP.) + ^[(B .V)B. + (B..V)Bo-V(Bo.B.)] 
+2p fK x n) + p ^ s - Po(^ • « s )^ - Po^f v s . 
We assume that (1) the accretion disc is infinitesimally thin, (2) the magnetic moment 
and the spin of the NS are parallel to the z axis, and normal to the disc, i.e., Bo = (0, 0, Bq) 
and f2 = (0, 0, Q) close to the inner edge of the disc, and (3) the x and y axes are along 
the disc plane, and the MHD wave is assumed to propagate in the xoz plane, i.e., the wave 
vector k = (ksin9, 0, kcos9), where 9 is the angle between the z axis and k. After carrying 
out Fourier transformation for Eqs. (14)-(16) we get the following dispersion equations, 
Q 2 ul k 2 V 2 A k 2 c 2 sin 2 9^ k 2 c 2 sin9cos9 2fli 

(1 + — 2 1 T- — ) V >* = —2 V " + V °V, 17 

„ n 2 u 2 k 2 Vjcos 2 9^ 2ili 

! + — --! — 2 K = V SX , (18) 



± 

UJ^ L) 2 OJ* UJ 

ujI k 2 c 2 cos 2 9, k 2 c 2 sin9cos9 



1 " —2 — >sz = ~^—2 Vs*, 19 

UJ 2 UJ 1 UJ 1 

where i is imaginary unit, V A (= yfB 2 /pp ), c s (= y/jPo/po), uj k (= yfGM/rf), and uj are 
Alfven velocity, acoustic velocity, Keplerian angular velocity, and angular velocity of the 
perturbation at r respectively, and v sx , v sy , v sz are the three components of the perturbed 
quantity of the speed. Equations (17)-(19) show that there exist three resonance MHD 
modes. At the magnetospheric radius r , the magnetic energy density is equal to the total 
kinetic energy density, i.e. B 2 /8ir = pV^/2 ~ pV k 2 /2 (Davidson & Ostriker 1973; Ghosh et 
al. 1977). Since the characteristic wavelength is in the same order with the magnetospheric 
radius (Rezania & Samson 2005), we then have JzVa ~ kVx ~ Vk/tq = uj k , or KVa = r)oj k . 
Because the thermal pressure of the plasma might be comparable with the magnetic pressure 
(c s ~ Va) just inside the magnetosphere (Miller et al. 1998), we also suppose kc s = \uj k . 
Here both rj and A are taken to be constant for certain sources. Substitute these relations 
into Eqs. (17)-(19) we can get the resonant modes. Specifically when 9 = 0, from Eqs. (17) 
and (18) we can get v sx = ±iv sy , i.e. v sx e lk ' v ~ lult = v sy e % ' r ~ lu) l z. Substituting this relation 
into the Eqs. (17)-(19) can give 



W^yJl + tfWk-n, (20) 

uj 2 = Vl + X 2 uJ k , (21) 
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LU 3 = ^Jl+T] 2 UJ k + VI, (22) 

where the negative solutions are excluded. Similarly when 9 = n/2, 

u{ = u>l (23) 
ul = <J* + fi 2 + ^( v 2 + A 2 ) + i| v /^V + A 2 ) 2 + 8fi 2 (^ + A 2 + 2), (24) 

ool = ujI + Q 2 + ^( V 2 + A 2 ) - ^uj 2 (^ + X^ + 8Q 2 ( V 2 + X 2 + 2). (25) 

Note that in the latter case the MHD wave couldn't propagate very far in the accretion 
disc, so we consider the coupling modes that propagate along the z axis (6 = 0) as a more 
promising explanation of the QPOs. 

We first rule out the possibility of the uj 2 mode as the lower kHz QPOs, since Eq. (21) 
requires that it should always higher than VL for stable accretion, which is contradicted 
with observations. The coupling between oj\ and u> 3 can also be excluded, which implies a 
constant frequency separation. So we suggest the upper and the lower kHz QPO frequencies 
be v-i = U2/2n and v\ = ui/2n. From Eqs. (20) and (21) we can get the following relation 
between the frequencies of the upper and the lower kHz QPOs, 



1 + A 2 , 

^ = )lTT7 {ui + Us) (26) 

where v s = Q/2n, or 

^ = ^=L=(^ + l)(whenr ? >A), (27) 

V -l = VTT5*(- + 1) (when v < A), (28) 

where e 2 = {jf - A 2 )/(l + A 2 ) and 5 2 = (\ 2 - rf) / {I + rf). Equations (27) and (28) indicate 
that the twin kHz QPOs may be divided into two groups, with the slope of the vz/vs vs- 
v\/v s relation either larger or smaller than 1. In the following we call them the large slope 
coefficient sources (LSCS) and the small slope coefficient sources (SSCS), respectively. 



3 COMPARISON WITH OBSERVATIONS 

We compare in Fig. 2 the z^/Vs vs - v \jv s relations obtained in last section with the observed 
kHz QPOs in six sources 4U 0614+09, 4U 1608-52, 4U 1636-53, 4U 1728-34, 4U 1915-05, 
XTE 1807—294, in which both the spin and twin kHz QPO frequencies have been measured. 
The spin frequencies, disposed in Table 1, are from van der Klis (2006), Mendez & Belloni 
(2007), Yin et al. (2007), Altamirano et al. (2008), and their references. For 4U 0614+09 we 
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adopt the updated spin frequency of 415 Hz (Strohmayer, Markwardt & Kuulkers 2008). The 
dots with error bars represent the measured values, and the solid lines stand for theoretical 
relations. We distinguish the vijv s vs. v\jv s relations for SSCS and LSCS, and accordingly 
adopt relation (27) to fit the data for 4U 0614+09, 4U 1608-52, 4U 1636-53, and 4U 
1728-34, and relation (28) for the other two sources, 4U 1915-05 and XTE 1807-294. For 
each source, the value of e or 8 for best fitting is also shown in the figure. It is noted that a 
cluster of the values (~ 0.3 — 0.9) of e and 8 can well reproduce the observed relations. In the 
left panel of Fig. 3 we show the observed and predicted relations for all of the six sources. 
In the right panel we plot the relation between Avjv s and v\jv s by use of the parameter e 
or 8 that we have got. 

For SSCS the peak separation of the twin kHz QPOs is less than the spin frequency, 
i.e., Av — v s = —(1 — l/Vl + s 2 )(vi + v s ) < 0, and decreases with z/i or z/ 2 ; for LSCS the 
peak separation is more than the spin frequency, i.e., Av = (v^l + 8 2 — + v s > u s , and 
increases with the increasing v\ or z/ 2 . In the former group, Av is around v s for 4U 0614+09 
and 4U 1728-34, and v s /2 for 4U 1608-52 and 4U 1636-53 (Miller et al. 1998; van der Klis 
1997; Stella, Vietri & Morsink 1999; Lewin & van der Klis 2006; M'endez & Belloni 2007). 

Our final note is that when the Alfven speed is equal to the acoustic speed of the plasma, 
i.e. r] = A, Eq. (26) will recover to the expression in the sonic-point beat-frequency model 
(Miller et al. 1998). In this case the peak separation is equal to the spin frequency and 
almost invariant. 

4 DISCUSSION AND CONCLUSIONS 

In this paper we propose a resonant MHD model for the twin kilohertz QPOs in LMXBs. 
The modes of the MHD waves vertical and parallel to the accretion disc are derived, and the 
twin kHz QPOs frequencies are identified with the frequencies of the two resonant modes. In 
this model the twin kHz QPO frequencies are correlated with the spin frequencies, and the 
separation frequencies also change with the QPO frequencies. We show that the measured 
relations between u 2 , and u s can be accounted for with reasonable values of the input 
parameters. 

There are several spin-involved MHD models for kHz QPOs in the literature. Osherovich 
& Titarchuk (1999) suggest that kHz QPOs can be explained as oscillations of large scale 
inhomogeneities (hot blobs) thrown into the NS magnetosphere. Participating in the radial 
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oscillations with Keplerian frequency, such blobs are simultaneously under the influence 
of the Coriolis force. The derived frequency relation is v\ = v\ + (2u s ) 2 , or (z/2/^) 2 = 
(v\lv s ) 2 + 4, which is plotted in the dotted curve in Fig. 4. The significant deviation from 
the measured data indicates that this model is not successful for most of the LMXBs. In 
their MHD loop oscillation model Li & Zhang (2005) derive a linear frequency relation, i.e., 
v< il v s — £,v\jv s -\- \ where £ ~ 1 is an input parameter. Here the upper kHz QPO frequency is 
assumed to the Keplerian frequency and the lower kHz QPO is identified as the principal fast 
kink mode of the standing MHD waves along the toroidal field lines at the magnetospheric 
radius. The relation is also plotted in Fig. 4 in solid curves. A comparison with the measured 
data shows that the fit is acceptable for all the six sources. 

Rezania & Samson (2005) propose a model for QPOs in LMXBs based on oscillating 
magnetohydrodynamic modes in NS magnetospheres. They argue that the interaction of the 
accretion disc with the magnetosphere can excite resonant shear Alfven waves in a region 
of enhanced density gradients, where accretion material flows along the magnetic field lines 
in the magnetosphere. The predicted v-ijv\ ratio is found to be independent of the NS spin 
frequencies. The main difference between Rezania & Samson (2005) and this work lies in that 
Rezania & Samson (2005) assume that the strong gravity of the NS produces a converging 
flow which will hit the star's magnetosphere in a large velocity, while we consider the motion 
of the plasma is still mainly Keplerian, before they enter the magnetosphere and corotate 
with the magnetosphere. We also include the effect of the gravity and the rotation of the 
NS for the trapped plasma, and find that the resulting wave frequencies relate with both 
the spin frequency and the gravity. 

So far we have assumed that the spin and magnetic axes of the NS are aligned in LMXBs. 
However, the existence of persistent millisecond pulsations in some LMXBs indicates that 
the magnetic axis is at least somewhat tilted from the spin axis. In this case the magnetic 
field is no longer homogeneous at a given radius in the accretion disc, and the orbit of plasma 
in the inner region of the disc becomes non circular. Furthermore, the stellar magnetic field 
can induce disk warping and precession (e.g. Lai 1999, 2003), and modulate the orbit and 
hence the QPO frequencies. For X-ray binaries the disc precession timescale is usually of 
tens of days to years (cf. Wijers & Pringle 1999 and references therein), which is much longer 
than the duration of each QPO observation. Additionally for accretion-powered millisecond 
pulsars the magnetic inclination is likely to be very small (Lamb et al., 2008). For the above 
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reasons we expect that the change of the QPO frequencies induced by oblique magnetic 
fields might be very small compared with the uncertainties in the measured frequencies. 

Our results indicate that the peak separation is always related to the spin frequency. 
What's more, there seems to be a weak positive correlation between the spin frequency 
and the parameter e for SSCS, which can be described as e = 2.28(±0.16)(z/ s /1000 Hz) — 
0.53(±0. 08), plotted in the left panel of Fig. 5. Substitute this relation into Av = (l/Vl + e 2 )(vi + 
v s ) — v\ we get a trend of Av changing with the spin frequency, as plotted in the right panel 
of Fig. 5. From the dark black curve to the light gray curve v\ changes from 1000 Hz to 100 
Hz in a step of 100 Hz. We find that when v s increases, Av varies from ~ v s /2 to ~ v s , and 
finally to ~ v s /2. The transitions occur at v s ~ 100 Hz and 500 Hz, respectively. 

The accretion process can take place only when the magnetospheric radius is less than 
the corotation radius (e.g. Ghosh & Lamb 1979). In other words, the Keplerian frequency 
at the magnetosphere radius should be more than the spin frequency of the NS if accretion 
process can take place. Because of this there is a minimum value of the lower frequency of 
the twin kHz QPOs in SSCS, v 1 > (y/1 + e 2 — l)v s . Besides, due to the fact that the peak 
separation must be positive values we can get the maximal value for the upper frequency, 
^2 < v s / (y/1 + e 2 — 1) for SSCS. These may serve as possible evidence to testify this model 
with future measurements of kHz QPOs in LMXBs. 
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Figure 1. The coordinate system centered at the magnetospheric radius. The x axis is along the radial direction, y axis the 
toroidal direction, and the z axis is normal to the accretion disc. 
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Figure 2. The relations between V2/v s and v\jv a for the six sources (for the measured data of 4U 0614+09: van Straaten 
ct al. 2000; van Straaten ct al. 2002; 4U 1608-52: van Straaten, van der Klis & Mendez 2003; 4U 1636-53: Altamirano et al., 
2008; Di Salvo, Mendez et & van der Klis 2003; Jonker, Mendez & van der Klis 2002; Wijands ct al. 1997; 4U 1728-34: Migliari, 
van der Klis & Fender 2003; Di Salvo ct al. 2001; Jonker, Mendez & van der Klis 2000; 4U 1915-05: Boirin et al. 2000; XTE 
1807-294: Linares et al. 2005; Zhang ct al. 2006). 
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Figure 3. Left The relation between vijvs and v\jv a . Right The relation between Au/u B and v\jv s . 



orange dotted curve Osherovich & Titarchuk; other lines Li & Zhang 
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Figure 4. The relations between V2/v a and v\/v s compared with the predicted ones in Li & Zhang (2005) and Osherovich & 
Titarchuk (1999). 




Figure 5. Left The relation between the parameter e and the NS spin frequencies in the four SSCS. Right The relation 
between the peak separation of the twin kHz QPOs and the spin frequencies with different lower kHz QPO frequencies. 
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Table 1. The measured and fitted parameters for six LMXBs. For the first four sources we use the relation U2 = (u\ + 
^s)/(Vl + e 2 )i and for the last two sources with V2 = Vl + S 2 (ui + u s ). 



sources 


v s (Hz) 


e or S 


error(±) 


X 2 /DoF 


slope 


4U f728-34 


363 


0.27386 


0.04629 


550.3/33 


0.96449 


4U 4608-52 


649 


0.89454 


0.00826 


450.5/16 


0.74531 


4U 4636-53 


584 


0.77732 


0.00768 


551.3/24 


0.78953 


4U 0644+09 


445 


0.45296 


0.04448 


563.0/39 


0.91091 



4U 4945-05 270 0.41657 0.01784 840.3/ 4 1.08330 

XTE 1807-294 190.6 0.33026 0.04243 94.6/12 1.05312 



